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Abstract 

Proper  orthogonal  decomposition  (POD)  technique  (or  the  Karhunan  Loeve  procedure) 
has  been  used  to  obtain  low  dimensional  dynamical  models  of  many  applications  in  en¬ 
gineering  and  science.  In  principle,  the  idea  is  to  start  with  an  ensemble  of  data,  called 
snapshots,  collected  from  an  experiment  or  a  numerical  procedure  of  a  physical  system.  The 
POD  technique  is  then  used  to  produce  a  set  of  basis  functions  which  spans  the  snapshot 
collection.  When  these  basis  functions  are  used  in  a  Galerkin  procedure,  they  yield  a  fi¬ 
nite  dimensional  dynamical  system  with  the  smallest  possible  degrees  of  freedom.  In  this 
context,  it  is  assumed  that  the  physical  system  has  a  mathematical  model,  which  may  not 
be  available  for  many  physical  and/or  industrial  applications.  In  this  paper,  we  consider 
the  steady  state  Rayleigh-Benard  convection  whose  mathematical  model  is  assumed  to  be 
unknown,  but  numerical  data  are  available.  The  aim  of  the  paper  is  to  show  that,  using  the 
obtained  ensemble  of  data,  POD  can  be  used  to  model  accurately  the  natural  convection. 
Furthermore,  this  approach  is  very  efficient  in  the  sense  that  it  uses  the  smallest  possible 
number  of  parameters  and  thus  is  suited  for  process  control.  Particularly,  we  consider  two 
boundary  control  problems:  (a)  tracking  problem  and  (b)  avoiding  hot  spot  in  a  certain 
region  of  the  domain. 
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Modeling  and  Control  of  Physical  Processes  using  POD 


1  Introduction 

There  are  many  families  of  patterns  for  which  it  is  possible  to  obtain  a  useful  systematic  char¬ 
acterization.  Often,  the  motivation  is  that  the  family  is  of  low  dimension,  that  is,  in  some  sense, 
any  given  member  of  the  family  might  be  represented  by  a  small  number  of  parameters.  These 
types  of  families  of  patterns  occur  frequently  in  both  nature  and  in  the  literature.  Such  examples 
include  turbulent  flows  [1,  2,  3,  4],  image  processing  [5,  6],  data  compression  [7],  human  speech  [8], 
and  human  faces  [9].  The  technique  applied  in  these  applications  is  known  as  the  proper  orthog¬ 
onal  decomposition  (POD),  although  it  also  goes  by  other  names  such  as  the  Karhunen-Loeve 
decomposition  [10],  principal  component  analysis  [11],  and  the  Hotelling  transform  [12]. 

POD  seems  to  have  been  first  proposed  [13]  in  1901  and  then  again  [14]  in  1933.  Since  its 
introduction,  the  method  has  received  much  attention  as  a  tool  to  analyze  complex  physical 
systems.  One  such  important  application  was  the  attraction  of  spatial  organized  motions  in  fluid 
flows.  Theodorsen  [15]  and  later  Townsend  [16]  observed  and  indicated  that  there  are  large-scale 
organized  motions  embedded  in  turbulent  shear  flows.  Lurnley  [1]  used  POD  as  an  unbiased  way 
for  extracting  structures  from  turbulent  flows.  Other  applications  of  POD  include  channel  flows 
[17,  18],  square-duct,  flows  [19],  and  shear  flows  [20,  21].  Other  scientists  have  also  applied  the 
POD  technique  to  fluid  related  problems.  For  instance,  POD  has  been  applied  to  the  Burgers’ 
equation  [22],  the  Ginzburg-Landau  equation  and  the  Benard  convection  [23].  Ly  and  Tran  [24] 
have  used  POD  to  simulate  and  solve  an  optimal  control  problem  for  fluid  flows  in  a  horizontal 
chemical  vapor  deposition  reactor.  In  addition,  Banks,  del  Rosario  and  Smith  [25]  have  applied 
POD  to  design  feedback  control  for  a  thin  shell  model.  The  goal  of  the  approach  is  to  represent 
an  ensemble  of  data  (called  snapshots),  obtained  from  physical  experiments  or  from  detailed 
numerical  simulations,  in  terms  of  an  optimal  coordinate  system.  That  is,  the  snapshots  can 
be  generated  by  a  smallest  possible  set  of  basis  functions.  Among  the  optimality  properties  is 
the  fact  that,  used  in  conjunction  with  the  Galerkin  method,  the  POD  yields  an  optimal  set  of 
basis  functions  in  the  sense  that  the  resulting  truncated  system  of  ODE’s  captures  the  maximum 
amount  of  kinetic  energy  among  all  possible  truncations  of  the  same  order. 

In  general,  to  implement  numerical  methods  such  as  POD  procedure  used  together  with  the 
Galerkin  method,  finite  element,  finite  volume,  or  spectral  method,  for  applications  in  engi¬ 
neering  or  science  it  is  assumed  that  the  physical  system  can  be  described  by  a  mathematical 
model  (in  general,  a  nonlinear  distributed  parameter  system  with  boundary/initial  conditions). 
The  governing  equations  are  usually  derived  based  on  some  simplified  assumptions  about  the 
physical  system.  On  the  other  hand,  many  applications  in  engineering  or  science  are  complex 
enough  that  a  mathematical  description  is  either  not  possible  or  when  a  model  is  obtained  it 
can’t  be  implemented  without  undue  complications.  For  example,  the  degree  of  freedom  of  dis¬ 
tributed  parameter  systems  is  essentially  infinite  and  the  relevant  mathematical  theory  is  either 
too  complicated  or  is  still  lacking  (in  the  case  of  nonlinear  distributed  parameter  system)  to  be 
implemented  in  an  industry.  In  this  study,  we  consider  a  physical  process,  whose  mathematical 
model  is  assumed  not  to  be  available,  but  data  can  be  collected  from  the  system.  For  example, 
in  a  physical  setting,  data  can  be  collected  from  thermal  images  or  from  thermocouples.  Our  pri¬ 
mary  goal  is  to  show  that  the  characteristic  features,  extracted  from  the  collected  data  using  the 
POD  procedure,  can  be  used  to  model  the  physical  system  efficiently  and  accurately.  Moreover, 
due  the  optimality  of  the  POD  technique,  this  approach  guarantees  that  the  number  of  parame¬ 
ters  used  to  represent  the  physical  process  is  smallest  possible.  Consequently,  process  control  for 
such  system  becomes  apparent.  For  our  physical  application,  we  consider  the  Ra.yleigh-Benard 
convection  problem.  Here,  the  governing  equations  are  used  only  for  the  purpose  of  generating 
an  ensemble  of  data.  The  classical  Rayleigh-Benard  problem  offers  a  first  approach  to  highly 
complicated  convection  flow.  Natural  convection  is  both  a  challenging  and  an  interesting  subject 
due  to  the  coupling  of  the  fluid  flow  and  the  energy  transport.  It  has  been  extensively  studied 
in  the  literature  because  of  its  importance  to  many  engineering  applications  and  because  it  is 
often  found  to  be  the  controlling  mechanism  in  many  natural  occuring  processes.  The  numerous 
applications  are  found,  for  instance,  in  meteorology,  geophysics,  crystal  growth  in  microprocessor 
manufacturing,  and  building  heat  transfer.  An  excellent  overview  of  this  subject  can  be  found 
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in  [26]. 

More  specifically,  we  consider  in  this  paper  a  mixed  convection  in  a  square  cavity  heated 
on  one  wall  and  cooled  on  the  opposite  wall.  It  is  well  known  that  in  such  a  configuration,  for 
moderate  Rayleigh  number,  the  fluid  flow  is  characterized  by  a  single  large  cell;  the  fluid  ascends 
along  the  hot  wall,  descends  along  the  cold  wall,  and  turns  at  the  opposite  end  of  the  cavity. 
For  small  Rayleigh  numbers  the  flow  is  weak  and  heat  is  transferred  by  conduction  across  the 
fluid.  However,  with  increasing  Rayleigh  number  this  flow  becomes  unstable  (see,  e.g.,  [27]). 
As  noted  above,  thermal  convection  is  the  result  of  the  coupling  between  the  fluid  flow  and 
heat  transport.  However,  in  this  study,  the  collected  ensemble  of  data  consists  of  only  partial 
information,  namely,  the  thermal  data.  Using  this  ensemble  of  thermal  numerical  data,  the  POD 
expansion  produces  an  “ optimal ”  compressed  description  of  that  ensemble.  This  description  is  a 
truncated  series  expansion  in  terms  of  a  set  of  basis  functions.  As  an  evaluation  of  the  success 
of  the  procedure,  we  project  data  from  outside  of  the  set  of  snapshots  onto  the  set  of  optimal 
basis  functions.  Even  in  this  case,  the  relative  errors  were  remarkably  good  (less  than  3%).  In 
addition,  we  also  show  that  the  resulting  eigenfunctions  and  coefficients  associated  with  each 
eigenfunction  can  be  used  to  solve  effectively  two  boundary  control  problems.  One  application  is 
the  tracking  problem  of  finding  the  thermal  boundary  condition  so  that  the  isotherms  inside  the 
cavity  match  with  a  prescribed  temperature  distribution.  The  other  problem  is  to  determine  the 
hot  wall  temperature  so  that  certain  region  inside  the  cavity  stayed  below  a  given  temperature. 
These  studies  offer  an  important  first  step  to  modeling  and  real-time  control  of  complex  physical 
systems. 

The  organization  of  the  paper  is  as  follow.  We  describe  in  §2  the  proper  orthogonal  decompo¬ 
sition  procedure.  The  physical  problem  and  snapshot  construction  are  described  in  §3.  Section 
4  describes  the  computation  of  the  POD  basis  functions,  their  uses  in  modeling  convective  mo¬ 
tion,  and  error  analysis.  Finally,  applications  to  two  practical  boundary  control  problems  are 
discussed  in  §5. 

2  The  Proper  Orthogonal  Decomposition 

Let  Ui(x),  i  =  1,2,...,  A  denote  the  set  of  N  observations  (also  called  snapshots )  of  some 
physical  process  taken  at  position  x  =  (x,y).  The  average  of  the  ensemble  of  snapshots  is  given 
by' 

1  N 

[/=<  1/  >=-][>;  (x)  (2.1) 

i  —  1 

We  form  new  ensemble  by  focusing  on  deviations  from  the  mean  as  follows: 

Vi  =Ui-U  (2.2) 

We  wish  to  find  an  optimal  compressed  description  of  the  sequence  of  data  (2.2).  One  description 
of  the  process  is  a  series  expansion  in  terms  of  a  set  of  basis  functions.  Intuitively,  the  basis  func¬ 
tions  should  in  some  sense  be  representative  of  the  members  of  the  ensemble.  Such  a  coordinate 
system,  which  possesses  several  optimality  properties  (to  be  discussed  in  the  sequel),  is  provided 
by  the  Karim  lien  l.oeve  expansion,  where  the  basis  functions  $  are,  in  fact,  admixtures  of  the 
snapshots  and  are  given  by: 

N 

*  =  X^-(*)-  (2.3) 

i—1 

Here,  the  coefficients  are  to  be  determined  so  that  $  given  by  (2.3)  will  most  resembles  the 
ensemble  {V^x)}^ .  More  specifically,  we  look  for  a  function  $  to  maximize 

i—1 


(2.4) 
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subject  to 

($,$)  =  ||$||2  =  1, 

where  (-,  •)  and  ||  •  ||  denote  the  usual  L2  inner  product  and  L2-norm,  respectively. 

It  follows  that  (see,  e.g.,  [24] )  the  basis  functions  are  the  eigenfunctions  of  the  integral  equation 


where  the  kernel  is  given  by 


N 


c(x,x)  =  —  ^(x)Vi(x  ) 


i  —  1 


Substituting  (2.3)  into  (2.5)  yields  the  eigenvalue  problem: 


N 


^  '  L/ j  Oj  —  A 0 j 
3  = 1 


(2.5) 

(2.6) 

(2.7) 


where  Ly  =  Vj)  is  a  symmetric  and  nonnegative  matrix.  Thus,  we  see  that  our  problem 

amounts  to  solving  for  the  eigenvectors  of  an  N  x  N  matrix  where  N  is  the  size  of  the  ensemble 
of  snapshots.  Straightforward  calculation  (see  also  [24])  shows  that  the  cost  functional 


l^|(Ti,$)|2  =  (A$,$)=A. 


which  is  maximized  when  the  coefficients  s  of  (2.3)  are  the  elements  of  the  eigenvector  corre¬ 
sponding  to  the  largest  eigenvalue  of  L. 

Remarks.  Alternatively,  we  can  consider  the  discrete  Karhunen-Loeve  expansion  to  find  an 
optimal  representation  of  the  ensemble  of  snapshots.  In  the  two-dimensional  case,  each  sample 
of  snapshots  Ui(x,y)  (defined  on  a  set  of  n  x  n  nodal  points  (x,y))  can  be  expressed  as  an  n2 
dimensional  vector  u,;  as  follows: 

Uil 

U*2 


(2.8) 


L  ^'n2  J 

where  ii,;j  denotes  the  jt  li  component  of  the  vector  u; .  Here,  the  discrete  covariance  matrix  of 
the  ensemble  u  is  defined  as  [12]: 


Cu  =  E{( u  -  m„)(u  -  m,i)7  }, 


(2.9) 


where 

=  Tju}  (2.10) 

is  the  mean  vector,  E  is  the  expected  value.  Equations  (2.9)  and  (2.10)  can  be  replaced  by: 


C'a  = 


N 


'  N 

E- 

Li=l 


and 


1 

N 


i  —  1 

respectively.  The  size  of  the  covariance  matrix  is  now  of  dimension  n 2  x  n2,  where  n2  is  the 
number  of  nodal  points  (or  pixels  of  the  samples).  The  eigenvector  corresponding  to  the  largest 
eigenvalue  of  C*a  is  the  dominant  POD  mode.  However,  this  approach  is  inhibitely  expensive, 
especially  when  the  samples  are  of  higher  dimensions,  because  one  has  to  solve  a  very  large 
eigenvalue  problem. 
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3  Rayleigh-Benard  Convection  and  Snapshots’  Construc¬ 
tion 

In  this  section,  we  consider  the  steady  state  problem  of  Rayleigh-Benard  system  in  a  square 
cavity.  It  is  well  known  that  if  the  Rayleigh  number  exceeds  a  certain  critical  value,  buoyancy- 
driven  thermal  cellular  convection,  combined  with  an  increase  in  the  heat  flux  through  the  fluid 
layer,  will  begin.  The  mathematical  formulation  described  below  is  used  solely  for  the  purpose 
of  generating  the  numerical  data  ( snapshots )  to  be  used  in  the  POD  procedure. 

3.1  Rayleigh-Benard  Convection 

Let  Q  be  the  rectangular  cavity  { (a? ,y)  |  0  <  x  <  L,  0  <  y  <  L}  in  TR?  which  is  filled  with 
fluid.  The  top,  rtop,  and  bottom,  l'h„,,  walls  are  insulated,  the  left  wall,  rielt,  is  maintained  at 
a  constant  cold  temperature  To  while  the  right  wall,  rrigll, ,  is  heated  with  a  parabolic  profile 
4(1  'i  -1'  ‘  !J(I‘  ~  V )  +T0,  where  Tx  >  T0  (see  Figure  1). 


u=v=0 


Figure  1:  Dimensionalized  Problem  Definition  of  Rayleigh-Benard  Convection 


We  assume  that  the  fluid  is  viscous,  incompressible  and  that  the  flow  is  laminar.  Under  these 
assumptions,  the  governing  equations  for  the  system  are  given  by  (see,  e.g.,  [26]): 


Continuity: 

o' 

II 

> 

(3.1) 

Momentum: 

P ^  +  pu  ■  Yu  =  -VP  +  pgjftT  -  To)  +  /iY2u, 

(3.2) 

Energy: 

<9T 

Pcp  +  Pcpu  '  v T  =  «V“T, 

(3.3) 

where  u  =  [«,  i>]T,  T  and  p  represent  the  fluid  velocity,  temperature  and  pressure  respectively; 
p  stands  for  the  dynamic  viscosity,  k  the  thermal  conductivity,  cp  the  specific  heat,  p  the  density, 
j3  the  volume  expansion  constant,  g  the  gravity  and  j  the  vertical  unit  vector. 

The  hydrodynamic  and  thermal  boundary  conditions  on  the  fixed  surfaces  are: 
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At  the  top  and  bottom  walls  of  the  cavity 


u  —  v  =  0 
dT 


(3-4) 


dn 


=  0, 


•  At  the  left  wall 


At  the  right  wall 


T  = 


u  =  v  —  0 
T  =  T0, 

u  =  v  —  0 

4(Ti  -  To) 


(3.5) 


(3-6) 


L2 


y(L  -  y)  +  To, 


The  initial  conditions  on  are: 


u(a-,  y,  0)  =  u°(a:,  y) 
T(x,  y,  0)  =  T°(x,y). 


(3.7) 


We  next  rescale  the  equations  (3.1)  (3.3)  to  obtain  the  dimensionless  forms.  Let’s  denote  by 


Ra  = 


cPp2pg{T j  -  Tq)L3 

jJK 


the  Rayleigh  number  and 


p  Vcp 
Jrr  =  — — 

K 

the  Prandtl  number.  We  now  substitute  the  following  relations 


_  K  V Ra.Pr  - 

u  =  - u 

pCp  L 

p  _  paVRaPr p 
pepL 2 

x  =  Lx 


T  =  {TS-  T0)T  +  T0 

t  =  PcpL~  j 
K  V  RaPr 
V  =  Ly 


ill  equations  (3.1)-(3.3)  and  in  the  boundary  and  initial  conditions  (3.4)-(3.7),  then  rename 
the  variables  [u,  T,  P,  x,  y,  t]  as  [u,  T,  P,  x,  y,  t]  respectively,  to  obtain  the  following  non- 
dimensionalized  system  of  governing  equations: 


V  •  u  =  0, 


(3.8) 


7— +  7u  •  Vu  = -VP  +  7j'T  +  V2u,  (3.9) 

dT 

7—  +  7u-  VT=  V2T,  (3.10) 

dT 

u  =  v=  — ■  =  0  on  r,op  and  Tbot,  (3.11) 

On 

u  =  v  =  0,  T  =  0  on  rieft,  (3.12) 

u  =  v  =  0,  T  =  4y(l  -  y)  on  riigM,,  (3.13) 
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ii(a-,  y,  0)  = 


T(x,  y,  0)  = 


Pcp  L  o 


RaPr 


u  (x,y) 


1 


(3-14) 


Ti  -  To 


(TV,*/) -To), 


where  7  =  y  "pp-  The  domain  now  becomes  a  unit  square  cavity.  Consequently,  the  steady-state 
system  of  governing  equations  are  given  by  (see  Figure  2): 


V  -  u  =  0, 

711  •  Vu  =  —VP  +  7  jT  +  V2u, 

711  •  VT  =  V  2T, 
dT 

u  =  v  ~  =  0  on  rtop  and  rbot, 

On 

u  =  v  —  0,  T  =  0  on  r 
"  •  -  0.  T  =  4t/(  1  —  y)  on  riigM. 


(3-15) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 


u=v=0 


Figure  2:  Non-dimensionalized  Problem  Definition  of  Rayleigh-Benard  Convection 


3.2  Construction  of  Snapshots 

We  assume  that  the  temperature  on  the  left-hand-side  of  the  cavity  is  cold  at  300° K  and  the 
fluid  has  the  following  properties: 

<*=  5.24  x  10 ^97x10-^ 

P  1-0£t  k  =  580  x  10-3^ 

(3  =  1010  x  10“6T  L  =  0.10 ra, 

Then  for  each  heat  intensity  T\  applied  to  the  right  wall  of  the  cavity,  there  are  two  corresponding 
dimensionless  numbers,  the  Prandtl  number  Pr  and  the  Rayleigh  number  Ra.  And  for  the  above 
parameters,  these  constants  have  values: 

Pr=P£r_  =  q  876 

K 
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and 

Ra  -  CpP  ‘8<jL  m  -  300)  =  921.9(Ti  -  300). 

/IK 

It  is  well  known  [26]  that  when  the  Rayleigh  number  Ra  (and  equivalently  Tf)  exceeds  a  certain 
critical  value,  buoyancy  driven  thermal  convection  cells  will  begin  to  form  in  the  cavity.  In  this 
paper,  we  restrict  our  flow  to  laminar  regime  and,  consequently,  we  restrict  Rayleigh  numbers  to 
be  less  than  106.  This  corresponds  to  values  of  the  right  wall  temperature  T)  in  the  range  from 
300°  K  to  1385° K. 

The  ensemble  of  snapshots  (see  Figure  3)  used  in  this  paper  was  obtained  numerically  by 
utilizing  a  commercial  fluid  dynamics  package  called  FIDAP  version  7.6.  The  steady-state  system 
of  non-dimensionalized  equations  (3.15)-(3.20)  was  discretized  using  81  quadrilateral  (9-nodal 
quadratic)  elements.  The  resulting  discrete  nonlinear  equations  were  solved  for  the  velocity 
fields,  temperature  distributions,  and  pressure  using  values  of  Rayleigh  numbers,  Ra ,  equal  to 
1,000,  5,000,  and  from  10,  000  ..  .200,  000  in  increments  of  10,000.  The  corresponding  values 
of  the  right  wall  temperature,  i|,  ranges  from  301.08°A'  to  516.94°A.  Only  22  snapshots  of 
temperature  distributions  were  used  in  the  POD  procedure.  In  addition,  we  generate  four  extra 
temperature  snapshots  corresponding  to  Ra  equals  to  500,  71,250,  115,000,  and  250,000  for 
reasons  to  be  discussed  later. 


Ra=1.0E3 


Ra=6.0E4 


Ra=1.3E5 


Ra=2.0E5 


Figure  3:  A  sample  of  Calculated  Snapshots  of  Temperature  Distribution 


4  Construction  of  POD  Basis  Functions 

Let  [//^(x),  i  =  1,2,..,,  22  denote  the  set  of  22  temperature  snapshots  as  described  in  §3.2 
above.  Here,  we  use  the  superscript  Ra  to  denote  the  dependency  of  the  snapshots  on  the 
Rayleigh  number.  For  convenience  in  notation,  the  superscript  is  sometimes  omitted  in  the 
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following  discussion.  The  average  of  the  ensemble  of  snapshots  is  given  by 

1  22 

U=<U>=-J2utRa(x)-  (4-1) 

We  form  new  ensemble  by  focusing  on  deviations  from  the  mean  as  follows: 

Vi  -  Ui  -  U,  i  =  1,2,... ,22.  (4.2) 

We  then  apply  the  POD  procedure  to  fluctuations  of  the  temperature  from  the  mean,  Vi’s.  The 
coefficients  a.i  s  of  the  POD  basis  functions  f&j,,  defined  by  (2.3),  are  the  entries  of  the  fcth 
eigenvectors  corresponding  to  the  Mh  largest  eigenvalue  of  the  N  x  N  covariance  matrix  L  (see 
§2): 

Lij  =  ^2  In  Vi^Vj^  d*' 

Since  the  covariance  matrix  is  a  symmetric  and  nonnegative  matrix,  the  eigenvalues,  A*  are  real 
and  nonnegative.  We  arrange  the  eigenvalues  in  decreasing  order  as  Ai  >  At  >  •  •  •  >  A22  >  0. 
Thus  $1  is  the  basis  function  corresponding  to  the  largest  eigenvalue.  A  representative  of  these 
basis  functions  is  given  in  Figure  4. 


POD  Basis  #1 


POD  Basis  #2 


POD  Basis  #3 


POD  Basis  #4 


Figure  4:  A  sample  of  Calculated  POD  Basis  Functions 


Any  snapshot  in  the  ensemble  can  be  represented  as  a  linear  combination  of  the  basis  functions 
as  follows: 

N 

URa  =U  +YJ<*£a$k 

k  =  1 


(4.3) 
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where 

o£a  =  (9k,URa -U).  (4.4) 

and  N  =  ‘22. 

We  next  investigate  at  how  much  error  is  introduced  by  truncating  the  series  (4.3).  That  is, 
we  now  consider  the  approximation 


U 


Ra 

M 


U  ■ 


M 

E 

k  =  1 


_  Ra  /f, 
OLu  <3>k 


(4.5) 


where  M  <C  N (=  22).  It  is  well  known  that  (see,  e.g.,  [24])  the  expansion  (4.5)  is  optimal.  In 
particular,  among  all  linear  combinations,  the  POD  is  the  most  efficient,  in  the  sense  that,  for  a 
given  number  of  modes  M ,  the  POD  decomposition  will  capture  the  most  possible  kinetic  energy. 
In  fact,  one  way  to  quantify  the  error  of  the  approximation  is  by  computing  the  relative  error  as: 


Em 


II  u-u$ia 

IK'll 


(4.6) 


Figure  5  shows  the  relative  errors  of  all  22  snapshots  against  reconstructed  solutions  using 
M  POD  modes  where  M  =  1,  2,  3,  and  4.  We  note  that  for  M  =  3  the  relative  error  is  already 
less  than  0.5%. 


Figure  5: 


Relative  Errors  of  Each  Snapshots  Plotted  for  Various  Values  of  M 


The  crucial  question  is,  therefore,  how  to  choose  M ,  a  priori,  so  that  (4.5)  is  a  good  approx¬ 
imation  to  (4.3).  As  already  indicated  in  several  publications  (see,  e.g.,  [24]),  where 

\  denotes  the  vt  li  eigenvalue  of  the  covariance  matrix  L,  represents  the  average  kinetic  energy 
contained  in  the  first  M  modes.  Hence,  to  capture  most  of  the  energy  of  the  system  contained 
in  the  N  POD  elements,  it  suffices  to  choose  M  so  that  A,;  1 11  fact,  the  ra¬ 

tio  Y?,  Xj/  Y?_,  A i  yields  the  percentage  of  the  total  kinetic  energy  in  the  N  POD  elements 
that  is  contained  in  the  first  M  POD  elements.  Since  the  eigenvalues  are  in  descending  order, 
Ai  >  A2  >  •  •  •  >  Aj\r  >  0,  one  can  reasonably  expect  to  achieve  a  high  percentage  of  the  total  ki¬ 
netic  energy  with  M  sufficiently  smaller  than  N.  In  our  case,  for  an  ensemble  of  22  snapshots,  an 
approximation  with  M  =  1  yielded  a  ratio  of  .9334.  However,  when  M  =  3,  this  ratio  increases 
to  0.9996  resulting  in  a  truly  significant  computational  efficiency. 

To  be  useful  in  developing  both  open  loop  and  feedback  control  strategies,  we  now  demon¬ 
strate  that  by  using  the  same  POD  basis  functions,  we  are  able  to  model  the  Rayleigh-Benard 
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system  for  Rayleigh  number  equals  to  500,  71,250,  and  250,000.  We  note  that  the  values  500  and 
250,000  are  not  in  the  range  used  orginally  to  compute  the  ensemble  of  22  snapshots.  To  com¬ 
pute  the  POD  approximation  for  each  new  Rayleigh  number,  we  fit,  for  each  mode  k,  the  original 
Rayleigh  dependent  coefficients  cxfa,  by  a  cubic  spline  interpolation.  We  then  evaluate  the  cubic 
spline  polynomial  at  Ra  =  500,  71,250,  and  250,000  to  obtain  the  corresponding  coefficients 
df-s.  The  approximations,  as  computed  by  the  POD  expansion  using  M  =  3,  are  compared 
against  the  FIDAP  solutions  in  Figure  6.  The  plots  show  a  remarkably  good  agreement.  Also 
Figure  7  shows  the  relative  errors  plotted  against  the  number  of  POD  modes  M .  For  M  =  3  the 
maximum  relative  error  is  already  less  than  3%.  The  approximation,  U  +  Ylt-i  is  then 

used  to  compute  the  POD  solutions. 


POD  SOLUTION 


FIDAP  SOLUTION 


FIDAP  SOLUTION 


FIDAP  SOLUTION 


Figure  6:  Comparison  Between  POD  and  FIDAP  Solutions  for  Rayleigh  Numbers  Not  Belonging 
to  the  Ensemble  of  Snapshots 


5  Boundary  Control  Applications 

Optimal  control  problems  have  been  a  subject  of  interest  to  experimenters  and  designers  for  a 
long  time.  More  recently,  flow  control  has  become  of  substantial  interests  to  mathematicians 
and  computational  scientists.  Perhaps  the  most  spectacular  example  of  successful  flow  control  is 
that  of  aerodynamic  control.  Here,  one  determines  wing  flaps,  throttle,  etc.,  so  that  an  aircraft 
executes  a  desired  maneuver.  In  this  and  numerous  other  examples,  the  control  is  computed 
using  either  sophisticated  mathematical  models  involving  partial  differential  equations  and/or 
ordinary  differential  equations  or  simple  empirical  rules.  In  this  section,  we  will  demonstrate 
the  effectiveness  of  using  the  POD  expansion  of  the  ensemble  of  data  to  control  the  temperature 
distribution  inside  the  cavity  without  the  utilization  of  a  mathematical  model. 
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Figure  7:  Relative  Errors  Plotted  Against  M  for  Rayleigh  Numbers  Not  Belonging  to  the  En¬ 
semble  of  Snapshots 


Desired  Ra 

Calculated  Ra  (using  3  POD  modes) 

Relative  Error 

500 

500 

1.504e-07 

71,500 

71,284 

4.835e-04 

115,000 

114,939 

5.317e-04 

250,000 

249,296 

2.817e-03 

Table  1:  Optimal  Solutions  to  the  Temperature  Tracking  Problem 

5.1  Tracking  Problem 

Let  /,;(x)  denote  a  desired  temperature  distribution  inside  the  cavity.  In  our  case,  Td  is  the 
numerical  solution  to  the  dimensionless  steady-state  system  (3.15)-(3.20)  using  FIDAP.  We  want 
to  find  the  right-wall  temperature  (or  equivalently  the  Rayleigh  number),  so  that  the  POD 
expansion  given  by 

3 

U3nU  +  J2ufa*C  (5-1) 

k  =  l 

is  “close”  to  Td .  This  boundary  control  problem  can  be  formulated  as  finding  the  Rayleigh 
number,  Ra ,  to  minimize  the  cost  functional 

J{Ra)=  [  \U3  —Td\2  dx, 

J  n 

where  U 3  is  given  by  equation  (5.1)  and  fl  is  the  cavity  domain.  This  unconstrained  optimization 
problem  is  solved  using  a  combination  of  golden  section  search  and  parabolic  interpolation. 
The  results  are  tabulated  in  table  1.  Here,  we  are  tracking  four  different  temperature  profiles 
corresponding  to  Rayleigh  numbers  equal  to  500,  71,500,  115,00  and  250,000. 

5.2  Avoiding  Hot  Spots  Problem 

In  many  applications  including  the  fabrication  of  microelectronic  circuits,  it  is  desirable  to  avoid 
“hot  spots”  along  bounding  surfaces,  i.e.,  places  where  the  temperature  peaks  might  occur  (see, 
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e.g.,  [28]).  This  is  necessary  in  order  to  avoid  melt-down  or  structural  failures.  In  our  case, 
we  want  to  find  the  right  wall  temperature  so  that  the  temperature  inside  the  shaded  region  in 
Figure  8  is  not  greater  than  the  dimensionless  temperature  0.1.  That  is,  we  find  the  temperature 
on  the  right  wall  (equivalently  the  Rayleigh  number)  so  that  to  minimize  the  objective 

Ns  R 

J(Ro)  =  ^  max{(/:((“(x,;)  -  0.1,0}+  y^, 

i  —  1 

where  {/^“(x,;)  are  the  values  of  the  temperature  of  the  computed  solution  using  3  POD  mode 
at  the  nodal  points  x, ,  which  is  inside  the  shaded  region  in  Figure  8.  Ns  —  120  is  the  number 
of  nodal  points  inside  the  shaded  region.  The  term  is  included  in  the  cost  to  achieve  the 
balance  between  the  desired  objective  and  the  expended  energy  to  obtain  the  goal.  The  results 
are  plotted  in  Figure  9.  The  optimal  value  of  the  Rayleigh  number  is  8046.6893  and  is  compared 
against  two  non-optimal  solutions.  We  note  that  in  the  non-optimal  solution  corresponding  to 
Ha  —  3056.6893  the  right  hand  corner  (darker  shade)  is  above  the  dimensionless  temperature 
0.1.  At  the  other  extreme,  namely  Ra  =  13046.6893,  the  dimensionless  temperature  is  below  0.1 
but  at  the  expense  of  a  higher  heat  input  at  the  right  wall. 


Figure  8:  Avoiding  Hot  Spots  Problem  Formulation 


14 


Modeling  and  Control  of  Physical  Processes  using  POD 


NON-OPTIMAL  SOLUTION 


OPTIMAL  SOLUTION 


NON-OPTIMAL  SOLUTION 


Figure  9:  Avoiding  Hot  Spots  Optimal  Solution  Compared  Against  Two  Non-Optimal  Solutions 
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